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gravity theory is considered. The equations of motion for R and H are solved ana- 
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is calculated in one-loop approximation and the back-reaction of particle production 
on the evolution of R is taken into account. Possible implications of the model for 
cosmological creation of non-thermal dark matter is discussed. 
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1 Introduction 



During the last two decades the accumulated astronomical data have unambiguously proved 
that the universe expansion is accelerated. The data include the observation of the large 
scale structure of the universe, the measurements of the angular fluctuations of the cos- 
mic microwave background radiation, the determination of the universe age (for a review 
see p-J), and especially the discovery of the dimming of distant Supernovae |2j. 

The driving force behind this accelerated expansion in unknown. Among possible 
explanations, the most popular is probably the assumption of a new (unknown) form of 
cosmological energy density with large negative pressure, P < —g/S, the so-called dark 
energy, for a review see e.g. [3]. 

A competing mechanism to describe the accelerated expansion is represented by gravity 
modifications at small curvature, the so-called /(i?)-gravity theories, as suggested in ref. 
In these theories the standard Einstein-Hilbert Lagrangian density, proportional to the 
scalar curvature R, is replaced by a function f{R), so the usual action of General Relativity 
acquires an additional term: 

^ = -T^ I d'xV^f{R) + S^ = -'^ [ d'xV^[R + F{R)] + S,n, (1) 

lOTT J lOTT J 

where mpi = 1.22 ■ 10^^ GeV is the Planck mass and Sm is the matter action. 

The original version of such models suffers from a strong instability in presence of 
gravitating bodies [5] and because of that more complicated functions F{R) have been 
proposed P-[9], which are free from the mentioned exponential instability. It was claimed 
in an earlier paper [10] that the instability could be eliminated by a mere addition of terms 
of the type R'^/mq^'^~^^ to the action, where is a constant parameter with dimension 
of mass. To some extent such terms may suppress the instability for systems with rel- 
atively high mass density, g > Ig/cm^, but for less dense systems e.g. for those with 
g ~ 10~^'*g/cm^ this mechanism demands too large coefficients in front of i?^ (small m^), 
which are at odds with big bang nucleosynthesis. 

Though free of instability [5j, the models proposed in ^6n8j possess another troublesome 
feature, namely in a cosmological situation they should evolve from a singular state in the 
past [llj. Moreover, it was found in refs. [T2l[T3] that in presence of matter, a singularity 
may arise in the future if the matter density rises with time; such future singularity is 
unavoidable, regardless of the initial conditions, and is reached in a time which is much 
shorter than the cosmological one. In the standard Friedmann cosmology, in which the 
energy density decreases with time, a future singularity may appear in several scenarios 
(phantom cosmology, quintessence models, ...), but not in models |6l-E], as shown in 
ref. |13]. This statement is in straight disagreement with earlier papers [11]. 

The subject of modified gravity includes about 10'^ papers at the present time and it is 
impossible to quote many relevant works. For a review one may look into references [T5l - [T7] . 

The aforementioned problems can be cured by adding to the action an i?^-term, which 
prevents from the singular behavior both in the past and in the future. In the present work 
we study the cosmological evolution of the Universe in a theory with only an additional 
R^ term in the action, neglecting other terms which have been introduced to generate the 
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accelerated expansion in the contemporary universe. The impact of such terms is negh- 
gible in the hmit of sufficiently large curvature, \R\ 3> l-Rol; where Rq is the cosmological 
curvature at the present time. 

In other words, we study here the cosmological evolution of the early and not so early 
universe in the model with the action: 

where m is a constant parameter with dimension of mass. 

Cosmological models with an action quadratic in the curvature tensors were pioneered 
in ref. [18]. Such higher-order terms appear as a result of radiative corrections to the 
usual Einstein-Hilbert action after taking the expectation value of the energy-momentum 
tensor of matter in a curved background. In such models, like for instance the Starobinsky 
model [TS], the universe may have experienced an exponential (inflationary) expansion 
without invoking phase transitions in the very early universe. This model has a graceful 
exit to matter-dominated stage which is induced by the new scalar degree of freedom, the 
scalaron (curvature scalar), which becomes a dynamical field in i?^-theory. 

The reheating process, due to gravitational particle production from scalaron (curvature 
scalar) oscillations, leads to a transition to a Friedmann-like universe. These features of the 
model are thoroughly discussed, for instance, in refs. Cosmological dynamics of 

fourth-order gravity were investigated in several works, see e.g. [23] and references therein. 

A somewhat similar study was performed in ref. [23] where a version of massive Brans- 
Dicke (BD) theory without kinetic term (i.e. with BD parameter w = 0) was considered. 
The Hubble parameter and curvature demonstrate oscillating behavior which resembles 
the one found in our paper (and earlier in many others), but quantitative features are very 
much different. 

Beside i?^-terms, terms containing the Ricci tensor squared R^^R^'^ are induced by 
radiative corrections as well, and with similar magnitude. However, the natural magnitude 
of such radiatively induced terms is quite small. The characteristic mass parameter, in fact, 
is of the order of the Planck mass in both cases, which makes this situation non-interesting 
for applications discussed below. On the other hand, R^ cosmology (without R^yR'^^) has 
been considered in the literature with much larger magnitude of R^ than the natural value 
from radiative corrections. The assumption of large R^ terms is made ad hoc to formulate 
a model which could, for instance, cure singularities. We follow the spirit of those works. 
However, it could be worthwhile to study the consequences of more complicate models 
with both R^ and R^uR^"" terms. It may be a subject for future investigation. 

The paper is organized as follows. In sec. [2] the cosmological equations modified by 
the presence of i?^-term are presented and solved both analytically and numerically in the 
case of a radiation-dominated (RD) universe. 

We study two equivalent sets of equations, for the time-time component of the (modi- 
fied) Einstein equations and for their trace. In sec. [3] particle production by the external 
oscillating gravitational field is studied. First, we derive the equation of motion for the evo- 
lution of R with the account of the back-reaction from particle production. This leads to 
an exponential damping of the oscillating part of i?, while the non-oscillating " Friedmann" 
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part remains practically undisturbed. The particle production influx into the cosmologi- 
cal plasma is estimated in the case of a massless, minimally-coupled scalar fleld. In the 
conclusions we discuss possible implications of the scenario, in particular for heavy super- 
symmetric dark matter. 



2 Friedmann-like Universe in Gravity 

We consider the theory with action ([2]) and use the conventions r/^jy = diag(+, — , — , — ), 
— ^/sT^jy + ■ ■ ■ and R^u = R'^^au- We also use the natural units h = c = = 1 and 
deflne the Planck mass as m?pi = Gj/. 

The modifled Einstein equations for this theory read 

Rfiu - -^QtivR - (^Rfiu - -^RQtJiv + QpLuT^^ - T^piPv^ R = -^^Tf^u , (3) 

where = g^'^'D^Vy is the covariant D'Alembert operator. We assume the Friedmann- 
Robertson- Walker metric with the interval given by 



dr 



2 



1 — kr^ 



+ r'^dd'^ + sin^ d d^p^ 



(4) 



In what follows we will neglect the three-dimensional space curvatur^, hence setting = 0. 
The curvature scalar R is expressed through the Hubble parameter H = d/a as 

R = -QH - UH"^ . (5) 

Therefore, the time-time component of eq. ([3]) reads 

where over-dots denote derivative with respect to physical time t. 
Taking the trace of eq. ([3]) yields 

/2 + 3//i? + mM i? + -^T;" ) = 0. (7) 

This equation is a sort of Klein- Gordon equation for a homogeneous scalar fleld (the 
"scalaron") of mass m, with a source term proportional to the trace of the energy- 
momentum tensor of matter. The General Relativity (GR) case may be recovered when 
m — 7- oo. In this limit we expect to obtain the usual algebraic relation between the curva- 
ture scalar and the trace of the energy-momentum tensor of matter: 

= -SttT^ (8) 



^This is a very good approximation, at least during the radiation-dominated epoch. 
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However, unlike the usual GR, in higher-order theories curvature and matter are related to 
each other through a differential equation, not simply algebraically. Therefore, the theory 
may approach GR as m — )■ oo in a non-trivial way or even not approach it at all. 

For a perfect fluid with relativistic equation of state P = g/3, the trace of the energy- 
momentum tensor of matter vanishes and R satisfies the homogeneous equation. The 
GR solution R = satisfies this equation, but if one assumes that neither R nor R 
vanish initially, the general solution for R will be an oscillating function with a decreasing 
amplitude. The decrease of the amplitude is induced by the cosmological expansion (the 
second term in eq. ([7])) and by particle production by the oscillating gravitational field 
R{t). The latter is not included in this equation and will be taken into account below in 
sec. [31 

It can be easily shown that the left-hand side of eq. ([3]) is covariantly conserved, which 
in turn implies the covariant conservation of the energy-momentum tensor of matter. The 
latter allows to write the evolution equation for the matter content, assuming it to be a 
perfect fluid with energy density g and pressure P: 

g + 3H{g + P) = 0. (9) 

As is well known, only two of equations (E]), O, and are independent. 

From eq. (Q it follows that relativistic matter, having equation of state P = g/3, 
satisfies 

+ AHgn = . (10) 

In what follows we will use either the set of eqs. and (ITUi) or the set and ([7]) as 
the basic equations. They are of course equivalent but their numerical treatment may be 
somewhat different. 

There is a possibility of gravitational particle production, which may non-trivially affect 
the solutions of the above equations. In first approximation, however, we neglect such 
contributions, which will be dealt with later on in the final part of this paper. 

2.1 Evolution equations in dimensionless form 

During most of its history (in terms of redshift, not time) the universe was dominated by 
relativistic matter, or in other words, it was radiation-dominated (RD). An exception to 
the RD regime is of course the inflationary stage with vacuum-like energy density and the 
universe heating epoch at the end of inflation when most probably the (non-relativistic) 
matter dominated (MD) regime was realized. After that, relativistic matter was dominant 
till redshifts of order of 10^, which corresponds to the moment of the matter-radiation 
equality. From the observations of the abundances of light elements we know for sure that 
at the time of big bang nucleosynthesis (BBN) the universe was quite precisely dominated 
by relativistic matter. If earlier in the course of the universe expansion and cooling down 
there were first order phase transitions in the primeval plasma, relatively short periods 
of vacuum-like matter dominance could have taken place. We also expect some (small) 
corrections to the RD regime due to the existence of particles in the plasma with masses 
comparable to the universe temperature, and because of the trace anomaly in the energy- 
momentum tensor of matter which leads to 7^ even for massless particles. The 
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universe might have even been in a practically pure MD regime after the post-inflationary 
RD stage. This regime could have been created by primordial black holes which evaporated 
early enough to bring the universe back to the normal RD epoch [25]. The last possibility 
is especially interesting in the case of i?^-inflation since the initial oscillations of R would 
be damped due to particle production at the end of inflation, and the GR solution could 
be restored. All these deviations from the pure RD regime would in turn induce deviations 
from the GR solution {R = 0) and give rise to oscillations of R even if they were initially 
absent. 

In what follows, we study the cosmological evolution in the i?^-theory assuming rather 
general initial conditions for R and H and dominance of relativistic matter. The MD 
regime will be studied elsewhere. 

It is convenient to rewrite the equations in terms of the dimensionless quantities 
T = Hot, h = H/Hq, r = R/Hq, y = g / {?irn?piHl) , and u = m/HQ,wheTe Hq is the 
value of the Hubble parameter at some initial time to- Thus the following two equivalent 
systems of equations are obtained: 

^ +3^^-2)1+ Y^-°' (11) 
y' + Ahy = 

and 

r" + 3hr' + = , , , 

12 

r + Qh' + 12h^ = . 

Here prime indicates derivative with respect to dimensionless time r. If we impose the 
"natural" relativistic initial conditions: 

TO = 1/2, 

"■0 ~ ^ ' 



1/0 = 1 



we flnd that there exists the exact solution h = l/(2r), y = l/(4r^), and r = 0, which 
corresponds to the usual general relativity, but this solution may deviate from GR because 
of deviations of the real expansion regime from the purely relativistic one. In principle, 
depending on the initial conditions, the solutions may oscillate around some value of hr, 
which may itself be different from 1/2, and in fact this is what we will flnd later on. 

Any non-negligible deviation from the GR solution may lead to observable effects, and 
hence to observational constraints on m, the only free parameter of the model. We will 
tackle these problems below both analytically and numerically. 

2.1.1 Approximate Analytical Solutions 

First we assume that the deviations from GR are small and expand h = 1/ (2r) + hi and 
y = l/(4r^) + yi, assuming that hi/h -C 1 and yi/y <S 1, and linearize the system of 
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equations. It is convenient to introduce a new unknown function zi = h'l, so we obtain 
three first-order linear differential equations with time- dependent coefficients: 



2r 



— co'^ ] hi + Tuj'^yi 



h'l = Zl 



(14) 



1 2 
■hi - -yi 



T 



We can find an approximate analytical solution of this system in the limit of large times, 
or UT ^ 1. In this limit we can treat the coefficients as approximately constant and find 
the eigenvalues and eigenfunctions of the system of differential equations. This method 
essentially consists in separating "fast" and "slow" variables. 
The characteristic polynomial is 



and the eigenvalues (for large ojt) are approximately 



(15) 



r 



^2,3 ~ — — ±iuj . 



At 



The general solutions of the system ([T^ is a linear combination of eigenvectors Vj : 



[hi,zi,yi] = ^CjVjexp 



dr' Aj(r') 



where Vi = [1, -3/r, 1/r], 1/2,3 = [1, -3/(4r) izw, -l/(5r/4 ± iwr^)] and 



exp 



dr'XUr' 



and exp 



rfr'A2.3(r' 



„±iu)T /^3/4 

~ e n ■ 



(16) 



(17) 



Since the solution must be real, one should take its real and imaginary part. 

Let us note that eigenvectors Vj are normalized to almost constant values, up to terms 
of order 1/r^. In principle coefficients Cj depend upon time but this dependence is quite 
weak, Cj ~ C^o + Cji/r"^ and asymptotically negligible. 

The correction to the GR solution corresponding to the first eigenvalue Ai quickly 
decreases, since ~ 1/r^ and y^i^ ~ l/^^? and so it can be asymptotically neglected. 
The solutions for hi corresponding to A2,3 oscillate and decrease more slowly than the GR 
solution, in fact 

(2,3) sin(a;r + yg) 



hr 



n3/4 



(19) 



while the solution for the energy density also oscillates but drops down faster than the GR 
one, y ~ r^^^^^. The complete asymptotic solution for h has the form: 



h{r) 



1 Ci sin(co'r -|- ^9) 



2r 



n3/4 



(20) 
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For sufficiently large r the second term would start to dominate and the linear approxi- 
mation would no longer hold. Below we will obtain approximate analytical solutions even 
in the non-linear regime, in the high-frequency limit. 

However prior to that it would be instructive to find the solution for the equivalent set 
of equations and ([7]) in the same approximation of small deviation from GR. Defining 
as above, dimensionless time, r = Hot, dimensionless Hubble parameter, h = H/Ho, and 
dimensionless curvature r = R/H^, we rewrite these equations in the form: 



r = -6h' - 12h^, 
r" + 3hr' + = , 



(21) 



Introducing the new function qi = r', we obtain the system of three first-order linear 
differential equations: 

3 



Its characteristic polynomial is 

P(A) 



zr 

r' = qi 

1 2 

/ii = — r hi 

^ 6 t 



A + -1 [X' + ^+u 



(22) 



r 



2r 



so the approximate eigenvalues for ut ^ 1 are 

2 



Ai 



^2,3 ~ — ±iuj . 



T ' At 

As above, the solutions of the system (J22l) are linear combinations of eigenvectors Vj: 



(23) 



[hi,r,qi] = '^CjVjexp j (It' \j{T') 



(24) 



where 



and 



Vi = [1, 0, 0] , V^2,3 



1 . 3 ■ 

-—, — 77 — ^ r> 1, =t lUJ 

6(5/(4r) iiw)' ' 4r 



exp 



dr'Ai(r') 



and exp 



dr'A2,3(r') 



^±iu)T /^3/4 



(25) 



So the oscillating solution for hi is the same as that of eq. f[T9|) . However, the slowly 
varying (non-oscillating) solution for hi decreases more slowly, namely as instead of 

1/t\ 

Probably this difference is related to the freedom in the zeroth order GR solution, with 
respect to the shift of time: 

(26) 



2(r + 5) 2t 2r2 ' 
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Such freedom tells us that the terms of order 1 /r^ in the first-order corrections are in some 
sense arbitrary. So the solution /ii ~ is spurious and should be disregarded. Anyhow 
the non-oscillating solutions for hi quickly disappear asymptotically and can be neglected. 

The solutions found above describe the oscillations of the Hubble parameter around 
the GR value l/(2r). Moreover, the amplitude of such oscillations decreases more slowly 
than 1/r, so at some stage the oscillations will become large and the condition hi <^ h will 
no longer be satisfied. After this stage is reached, the linear approximation is no longer 
valid and the method used above becomes inapplicable. 

However, we can still find the asymptotics of the exact nonlinear equations ([5]) and ([7]) 
at large times looking for solutions in the form: 

h{T) = A{t) + Bs{t) sin z/r + Bc{t) cos ut , (27a) 

r(r) = C(r) + Ds{t) sin i/r + Dc{t) cos ut . (27b) 

Here coefficients A, B, C, and D are slowly varying functions of time and u, assumed large, 
may in principle be different from u. As we will see below, this is indeed the case due to 
radiative corrections (scalaron mass renormalization) . In the approximation taken here we 
find u = u. 

We obtain approximate equations for these functions equating the coefficients in front 
of the slowly varying terms, and in front of sinwr and coswr. These equations are ap- 
proximate because we do not take into account higher-frequency terms which appear as a 
result of non-linearity, but the taken approximation happens to be quite accurate. Doing 
so, we find from equations fl2Tl) : 

A' + 2^2 + B^^ + Bl = -C/6 , (28a) 
5; - B^u + AABs = -Ds/6 , (28b) 
B'^ + B,u + AAB^ = -DJ6 , (28c) 

C" + 3 (^AC + ^BsD'^ + ^B,D'^ - ^uBsD, + ^^vB.D^j + w^C = , (28d) 

- 2uD'^ - u'^D, + 3A(D^ - uDc) + 3C"5, + u'^D, = , (28e) 
D'^ + 2z/D; - u'^D, + 3A{D'^ + vD^) + 3C"5, + u'^D^ = . (28f ) 

Assuming 

A=-, 5. = ^, i?e=-, C=4, D, = ^, De=-, (29) 

r r r r r 

and keeping only the dominant terms (lowest powers of 1/r) we obtain the solutions 

= 00^ , (30a) 

ds = Gubc , (30b) 

dc = —6ubs , (30c) 

hl + hl = 2a{2a - 1) , (30d) 

c = 18a(l - 2a) . (30e) 



It is interesting that equation ( ]30dp demands a > 1/2 in presence of oscillations. We will 



see from the numerical solution that this is indeed the case. 
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To summarize, the situation is the following: we initially assume hi/h ^ 1 and solve 
the linearized systems or (122|) . finding that the Hubble parameter oscillates around 
the value l/(2r) with amplitude growing with time as hi/h ~ t^^^] eventually, such oscilla- 
tions (and hence non-linear terms) would become dominant and the linear approximation 
becomes invalid. However, we can proceed further using a sort of truncated Fourier expan- 
sion which allows to take into account the non-linearity of the system in the limit ur ^ 1. 
As a result we have found that hi/h — )■ const. In other words, the amplitude of the oscil- 
lating part of h asymptotically behaves as 1/r, i.e. in the same way as the slowly- varying 
part of h. 

To be completely sure about these analytical results we have to check if the exact nu- 
merical solution of the system (fTTj) shows the same behavior. Still, the analytical estimates 
presented above are of great interest for the calculation of the evolution of R and H when 
particle production effects are taken into account. 

2.1.2 Numerical Solutions 

We integrate the system of equations (fTTl) starting at tq = 1/2, with the initial conditions 

/lo = 1 + 

h'^ = -2 + 6h'^ (31) 
yo = l + 6yo, 

where 6ho, dh'^ and 5?/o do not vanish simultaneously. As expected, numerical integration 
with the initial conditions given by eq. f|T3|) gives the usual GR solution h = l/(2r) within 
numerical precision, so we are interested in the more general case in which the initial 
conditions deviate from the GR values. 

As it was said earlier, systems f lTT]) and f lT^ are equivalent. However, for the numerical 
integration of these systems one has to specify initial values of different quantities. For 
the integration of system (ITTj) one has to fix ho, /iq, and i/q, while for the integration of 
( !T2|) the values of h^, vq, and Tq must be specified. The expression of one set of initial 
values through the equivalent values of another set can be found using the equations under 
scrutiny. Indeed, once ho, and i/q are chosen, h'^ is uniquely determined through the first 
equation in (fTTll . and consequently tq and Tq are specified as well, via (fT2l) . After all, both 
systems are equivalent to the same single third-order differential equation, whose solution 
is determined by three initial conditions. 

We have found that the numerical solutions of the system of equations ( ITTl 131]) are 
in very good agreement with the previous analytical estimates in the linear regime, i.e. 
for initial conditions fulfilling the requirements 6ho/ho ^ 1, Syo/yo ^ 1, and S> 1. In 
figures [U [21 and [3] we present the numerical results for the dimensionless Hubble parameter 
h determined from the system of equations (fTTj) . with u = 10 and initial conditions 

Sho = 10-^ , Sh'^ = , 5yo = (fig. HD (32a) 
Sho = 10"^ , Sh'o = , 5yo = (fig. E]) (32b) 
5ho = 10~2 , 5ho = , Syo = (fig. E]) (32c) 
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The function hr is found to oscillate around the central value h = 1/2 with amplitude 
hi ~ r~^/^. As the deviation from the ideal GR behaviour increases, the average value of 
hr also varies, and in general it is no longer equal to 1/2. A very good functional form for 
fitting the solutions is 



hir) ~ ^ — = h hose , 



(33) 



where the dimensionless parameters a and /3 are very slowly varying functions of time. This 
fit is shown, for instance, in figures H] and O where the numerical solution for dh^ = 0.5 is 
presented. A deviation from the analytical estimates of the linearized equations is to be 
expected in this situation, since the condition Sho/ho ^ 1 is not fulfilled and non-linear 
terms in eqs. ffTTl) are important. 

For the moment, let us concentrate on the case of small 6 ho, in which we can safely 
take a = 1. Qualitatively, one notices that parameter /3, evaluated at the same value of 
UT, increases roughly linearly with the initial displacement 6hQ. Moreover, when solving 
the system of equations flTTl) with 6hQ = and varying Sh^, we find again a (roughly) linear 
relation of the form /3 oc 



h T 




^osc 

0.00006 
0.00004 
0.00002 

-0.00002 
-0.00004 
-0.00006 



h T^l^ 




COT 



Figure 1: 
and uj = 
6.29 X 10 



Numerical solution of eqs. (ITT]) with 6hQ 



10-^ 5h'Q 



0, yo 



10. The best fit, for functional form (!33|) . is given by a ~ 1, /3 

-5 



hT 



h T^'^ 




COT 




COT 



Figure 2: Numerical solution of system (ITT]) with dh^ = 10 ^, dh'^ = 0, yo = 1, 
CO = 10. The best fit is a ~ 1, /3 ~ 6.28 x 10"^. 



In figures [6] and [7] we present the numerical results for the initial conditions 

Sho = 1.5 , Sh'o = 0, 5yo = 0, u = 100 . (34) 
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LOT 



Figure 3: Numerical solution of system (ITT]) with 6hQ = 10 ^, (5/iq = 0, z/q = 1, 
00 = 10. The best fit is given by a ~ 1, /3 ~ 6.26 x 10^^. 




LOT 



Figure 4: Numerical solution of system (fTT]) with 6hQ = 0.5, Sh'^ = 0, i/q = 1, 
u = 10. The best fit at large times is a ~ 1.14 (dotted line in the left panel), 
P ^ 0.20. Please also note, in the right panel, that the oscillating part of h 
does not decrease exactly as r~^/^. The apparent up-down asymmetry in hose 
is due to the fact that a is a function of time, and that we centered oscillations 
around its value at late times. 



a 




LOT 



200 400 600 800 1000 



Figure 5: Evolution of a with time. Initial conditions are those of figure HI 



Results are, at least qualitatively, in agreement with the analytical estimates made in the 
non-linear regime, for ut ^ 1. Evidently, the amplitudes of the oscillating terms of both 
h and r decrease faster than r~^/^ (linear regime), and rather close to r~^. Furthermore, 
the Hubble parameter does not oscillate around the GR value hr = 1/2, but around a 
larger value, as expected from eq. fl30d|) . In fact, for these initial conditions, the best fit 
at the considered final integration time is given by a ~ 1.246. The universe evolution 
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Figure 6: Numerical solution for the dimensionless Hubble parameter h with 
initial conditions ( IMl) . Oscillations in the left panel are almost of constant 
amplitude, whereas in the right panel they are clearly decreasing. 




Figure 7: Numerical solution for the dimensionless scalar curvature r with 
initial conditions ( l34l) . 

in R + gravity after inflation, without relativistic matter, was considered in ref. |18] . 
where it was stated that the non-oscillating part of h tends asymptotically to the GR 
matter-domination value /ir = 2/3, in contrast to our relativistic case. We plan to study 
what happens in cosmology with non-relativistic matter in the future. 
Our result that hose ~ t~'^^^ f l20|) in the linear regime agrees with what found in the last 
paper of [IH], namely that hose ~ a~^/^. The latter, however, is not in perfect agreement 
with our results in the non-linear regime, where we have a ~ 1.246 but hose ~ 1/t. This is 
probably due to the fact that we have not been able to reach the true asymptotic regime, 
but only a pre-asymptotic region. 

3 Particle Production and Back-Reaction 

The particle production rate by the oscillating gravitational field in R"^ gravity was con- 
sidered in ref. [T9||20]. where it was estimated as F ~ m^/mp^. In this section we present 
more rigorous calculations, which are essentially in agreement with ref. [20] . We will derive 
below a closed equation of motion for the cosmological evolution of R with the account of 
the back-reaction of particle production. To this end we consider a massless scalar field 
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minimally-coupled to gravity. Its action can be written as: 



1 



In spatially-flat FRW background ^ it leads to the equation of motion: 

+ 3i^0 - ^ A0 = . 



(35) 



(36) 



Field enters the equation of motion for R ([7]) via the trace of its energy-momentum 
tensor: 

Tj:{<P) = -g^^d,<pd,ct^^-{d<pf. 

It is convenient to introduce the conformally rescaled field, x = ci{t)4'j aiid the conformal 
time 7], such that adr] = dt. It terms of these quantities we can rewrite the equations of 
motion as: 



„ a', „ „ 1 
R" + 2-R' + m^a^R = Svr^— 



a 



rripia^ 



,/2 
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(Vx)' + - -(xx' + x'x) 



R = -Qa"/a^ , 
X"- AX+^«'^X = 



while action fl55|) takes the form: 



dr]d'x x' -(Vx) - 



2 O.^R 2 



6 



-X 



(37) 

(38) 
(39) 

(40) 



Here and above prime denotes derivative with respect to conformal time. 

Our aim is to derive a closed equation for R taking the average value of the x-dependent 
quantum operators in the r.h.s. of eq. ( 1371) over vacuum, in presence of an external classical 
gravitational field R. Our arguments essentially repeat those of ref . [26] , where the equation 
was derived in one-loop approximation. 

We quantize the free field x^^"* as usual: 



d^k 



(27r)3 2Efc 



-ik-x _|_ ^ik-x 



(41) 



where x^ = [rj, x), /c^ = {E^, k), and k^k^ = 0. The creation-annihilation operators satisfy 
the usual Bose commutation relations: 



dk, d\ 



Equation (!39|l has the formal solution 



(42) 



X{x) = x^^Kx) - - I d'yG{x,y) a\y)R{y)x{y) = X^^K^) + ^xix) 



(43) 
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where the massless Green's function is 

G{x. y) = ((a^o - i/o) - |x - y|) = 7^^(A^ - r) . (44) 

47r|x — y| 47rr 

We assume that the particle production effects shghtly perturb the free solution, so that 5x 
can be considered small and the Dyson-like series can be truncated at first order, yielding 

^{x)^X^'\x)-\j d'yG{x,y)a\y)R{y)x^'\y)^x^'\x)+x^'\x). (45) 

We can now calculate the vacuum expectation values of the various terms in the right-hand 
side of equation (1371) . keeping only first-order terms in x^^^- All terms containing only x^^^ 
and its derivatives have nothing to do with particle production and can be re-absorbed by 
a renormalization procedure into the parameters of the theory, so they are of little interest 
here. The other terms are calculated using formulas such as 

j d^y G{x, y) a\y)R{y)x{y) = j A G{x, y) [a'RdyX + dy{a'R)x] , (46a) 

dk e*"^ = 7i6{a) + iV (^^^ . (46b) 

We have collected all the terms, arriving to the expressions: 

1 r\ ,a\r^')R{r^') 







^,^.^^,,^,^^^j\„.(mmr, ,4Tb) 

967r^ y^,, 7]-!]' 
1 f\ ,ia'iv')Riv')y 



{XX + x'x) ^ -TFT / dr]' . (47c) 

487r2 J 7] -7]' 



Substituting these expressions into (137|) . we obtain a closed integro-differential equation 
for R, of which we will find an approximate analytical solution. We also plan to find the 
numerical solution of this equation but this is a much more complicated problem. Still for 
our purpose the approximate analytical solution is accurate enough. 

First of all, one has to note that despite having oscillating H and R, the scale factor 
a basically follows a power-law expansion, so it varies very little during many oscillation 
times u~'^. Thus, we expect that dt]/-)] ~ dt/t and that the dominant part in the integrals 
in ( H5l) is given by derivatives of R, since R' ~ uR + t~^R ~ uR, because ut ^ 1. 
The dominant contribution of particle production is therefore given by eq. ( 147bl) . and yields 



U^SHH^m'H^-^^ir i,f (^mWr ^_^^f,,m, (48, 

12ir m|,a*/^ ' i, - r,' 12ir mj., t-t' ^ ' 

The equation is naturally non-local in time since the impact of particle production de- 
pends upon all the history of the evolution of the system. The equation is linear in R, 
in contradiction with reference [21], where the r.h.s. of the equation is quadratic in R. 
The latter result is physically doubtful because if the sign of R changes, the effect of the 
particle production term would not be a damping of oscillations but their amplification. 
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3.1 Analytical Calculations 

We repeat the calculations of section 12.1.11 using expansion fl27|) and including the back- 
reaction effects in the form of equation fHS]) . The right-hand side of this equation can be 
written as: 



to t t' Jq "T Je 'T Jo 

rt-to (j 

= g J dT—+ 

/•t—to 2 

-|- g cos(mit) / dr — [Fc cos(mir) — Fg sin(mir)] + 

Je ^ 

dr — [Fc sin(mir) + Fs cos(mir)] , (49) 



where 



1 ml 



12n rn?pi ' 



(50a) 

Fc = Dc + 2miDs - mlDc , (50b) 
Fs = Ds- 2mibc - m\Ds . (50c) 

To avoid confusion we need to mention that r here is simply the integration variable and 
is different from the dimensionless time r of the previous sections. 

We have introduced here the new notation mi which is equal to m plus radiative 
corrections specified below, and which corresponds to u in eqs. f l27al 12 7b I) . The difference 
between m and mi is not essential under the integral but it should be taken into account 
in the l.h.s. of eq. (37), where we should take mi instead of m. 

Please note that the slowly-varying functions C, Dg and Dc inside the integrals are to 
be evaluated at (t — r), and a dot denotes derivative with respect to (physical) time t, not 
r. Because of the 1/r factor, the integral is logarithmically divergent, but this divergence 
can be absorbed into the renormalization of mass, m, and coupling, g. So we separate 
the integral into two parts: one where r goes from to some small parameter e, which 
determines the normalization point at which the physical mass and coupling are fixed, and 
another, taken from e to (t — to)? which gives corrections to the physical qualities due to 
interactions. More details can be found in ref. |26] . 

Equating the coefficients multiplying the slow varying terms, sinmit, and cosmit in 
the same way as it has been done in sec. I2.1.H see eqs. fl28|) . we obtain the same first three 
equations fl28a[ I28bl I28cp . where the effects of particle production do not directly appear, 
and the remaining three ones with the additional terms coming from eq. fBHj) . see also 
expansion f H^ . The latter equations become integro-differential but they can be reduced 
to differential equations in the case of fast oscillations. So the complete set of equations 
with the account of particle production has the form (for convenience we also include the 
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unchanged first three equations of set (128|) : 

A + 2A^ + Bl + Bl = -C/6 , (51a) 

B, - B^nii + AABs = -DJ6 , (51b) 

Be + Bsfui + AABc = -DJ6 , (51c) 

C + 3AC+ IBsDs + Ib,D, - ImiBsD, + ^miB^Ds + m'C^g ! ' dr-, (51d) 

I I I I J t 

Ds + (m^ - ml)D, - 2miZ)c + 3A(Z), - miD^) + 

■ cos(mr) + sin(mr) 
+ 3Ci3s ~ fif / ar , (51e) 



T 

,2 



+ {m^ - mi)Dc + 2miDs + 3A{Dc + miDs) + 

+ 3CB. g dr ^' " ^' . (51f) 



r 



In integrals fl51ep and fl51fj) . containing quickly oscillating functions, the effective value of 
r is about 1/m. Thus we can approximate F{t — r) ^ F(t) and take such factors out of 
the integrals. Let us analyze now for example equation f l51ep term by term. The analysis 
of eq. ( I51f[) is similar. In what follows we neglect D in comparison with m?D. 



The dominant term in eq. f l51ep . which is the coefficient multiplying D^, determines 
the renormalization of m: 

ft-to J 

'2 '2 '2 I / \ 

nil = m +gm J — cos mr. (52) 

The next subdominant term, which is the coefficient in front of D^, determines the decay 
rate of Dc- 

■ qm /"*^*o nqm , , 

D, = ^D,J —sinmr^^D,. (53) 

We skipped here the term gDc, which leads to higher order corrections to the production 
rate. Thus the decay rate is 

r« = -^ = ^. (54) 

Correspondingly the oscillating part of i? or if behaves as 

cosmit — )■ e^^"^ cosmit. (55) 

We will use this result in the next subsection in the calculation of the energy density influx 
of the produced particles into the primeval plasma. 

3.2 Gravitational Particle Production 

From equation f p9|) follows that the amplitude of gravitational production of two identical 
X particles with momenta pi and p2 in the flrst order in perturbation theory is given by 

di]d^x {pi,P2\xx\0) , (56) 
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where the final two-particle state is defined by 



The factor 1 / -\/2 is simply the correct normalization of the two-particle state due to the 
Bose statistics. Using eq. (HT]) . we find 



= a/2 e*^"^''i'''"^''2)^-*(pi+P2)-x _ (57) 

Here i?| = k^, and the function a^R has the form 

a^{rf)R{vi) = Dijf) sin(a)r7) , 

where D{ri) is a slowly-varying function of (conformal) time, Cj is the frequency conjugated 
to conformal time. Under these approximations, the amplitude (156|) becomes 

A(pi,p2) = ^ /" drid^xD{ri) {e'^"^ - e"^^^) g^^^fi+^f^)" e-^^P^+P^^''' . (58) 
6v2 J 

Taking Ep^ > and neglecting at this stage the variation of D with time, we obtain 

A(pi,p2) ^ -^Dir])i2nY5^'\p, + p^) 5(i?p, + - u) . (59) 

In order to find the particle production rate per unit comoving volume and unit conformal 
time, we need to integrate the modulus squared of this amplitude over the phase space, 
namely: 



n' 



d'^Pid^P2 \A{pi,p2] 



{27r)H Ep,Ep, VAt] 



5767r 



(61) 



where V and At] are the total volume and conformal time, which of course go to infinity, 
n is the number density of the produced particles, and a prime denotes derivative with 
respect to conformal time. Since the energy of the produced particles is equal to uj/2, we 
find for the rate of gravitational energy transformation into elementary particles: 

n^^p^ (62) 
^ 2 11527r ^ ' 

and so the rate of variation of the physical energy density of the produced x-particles is 
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Here, (i?^) is the square of the amphtude of the oscillations of R and we substituted 
a) = am. To obtain the the total rate of the gravitational energy transformation into 
elementary particles we should multiply the above result by the number of the produced 
particle species, N^fj, so the total rate of production of matter is qpp = N^ffQ^. 
Note that result (1631) coincides with that of ref. [5D], although performed in a slightly 
different cosmological regime. 

Now we can calculate the evolution of the cosmological energy density of matter, which 
is determined by the equation: 

g = -AHq + Qpp . (64) 

We assumed here that the produced matter is relativistic and so the first term in the r.h.s. 
describes the usual cosmological red-shift, while the second term is the particle source from 
the oscillations of R. Since q is not oscillating but a smoothly varying function of time, its 
red-shift is predominantly determined by the non-oscillating part of the Hubble parameter. 
He ~ a/2t, see eq. (!33|) . 

Parametrizing the oscillating part of the Hubble parameter as Hose — P cos mt/t, we 
find for the oscillating part of curvature: 

„ 6/3msinmt p , , , 

R ~ — - — e-^«* . (65) 

Here we took into account the exponential damping of R, which was for brevity omitted 
in the expression for H just above. 

Correspondingly the energy density of matter obeys the equation: 

This equation can be explicitly integrated as it is, but for a simple analytical estimate we 
will use the instant decay approximation. Namely we neglect the exponential damping 
term, when 2rpt < 1, and take a = ai = 1.25, according to the numerical estimate of 
sec. 12.1.21 For 2Tpt > 1 we completely ignore the second (source) term in eq. (1^ and 
take a = 0^2 = 1. This choice corresponds to the GR solution and we believe that it is 
realized when the oscillations disappear, as follows from the analytical estimates presented 
above. Thus at short times, 2rpt < 1, the energy density of matter would be: 

" "'"{t ) + 32^(2ai - l)t \ J ■ ^"'f 

For large times, i.e. 2Tpt > 1, equation ( l66i) becomes homogeneous and its solution 
is simply the relativistically red-shifted energy density with the initial value determined 
from eq. §7^ at t = l/i2Tpy. 



Q 



768nm% (2r«t)2- [i (2^^-^^)^-"^ + £^ (l - (2r«t,„)'"^-') 
where we parametrized the energy density of matter at the initial time tin as 



(6^ 



in 
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Parameter k is arbitrary, and depends upon the thermal history of the universe before tin- 
la particular, k = is possible and does not contradict our picture, since the equations of 
motion have non-trivial oscillating solutions even if ^ = 0. 

The first term in equation fl68p is the contribution of normal thermalized relativistic 
matter, while the second also describes relativistic matter, but this matter might not be 
thermalized, at least during some cosmological period. Depending upon parameters the 
relative magnitude of non-thermalized matter might vary from negligibly small up to being 
the dominant one. 



4 Discussion and Implications 

The characteristic decay time of the oscillating curvature is 



1 24mL ^/lO^GeVA^ 



2Tii m 



\ m J 



^R=7^ = -zf^^M^:—] seconds. (70) 



The contribution of the produced particles into the total cosmological energy density 
reaches its maximum value at approximately this time. The ratio of the energy den- 
sity of the newly produced energetic particles and that of those already existing in the 
plasma, according to eq. is: 

gh^ _ 8(3^Neff 1 - (2rfiti„)2-i-i ^^^^ 



gtherm /^(2ai - 1) {2T nUn)^'^^-'^ 

If we take tj„ ^ 1/m, then tinT^ ~ m? /m?pj^ <^ 1 and the effects of non-thermalized matter 
may be negligible. However, for sufficiently large /3 and possibly small k the non-thermal 
particles may play a significant role in the cosmological history. 

The infiux of energetic protons and antiprotons could have an impact on BBN. Thus 
this would either allow to obtain some bounds on m or even to improve the agreement 
between the theoretical predictions for BBN and the measurements of primordial light 
nuclei abundances. 

The oscillating curvature might also be a source of dark matter in the form of heavy 
supersymmetric (SUSY) particles. Since the expected light SUSY particles have not yet 
been discovered at LHC, to some people supersymmetry somewhat lost its attractiveness. 
The contribution of the stable lightest SUSY particle into the cosmological energy is pro- 
portional to 

and for msusY in the range 100 — 1000 GeV the cosmological fraction of these particles 
would be of order unity. It is exactly what is necessary for dark matter. However, it 
excludes thermally produced LSP's if they are much heavier. If LSP's came from the 
decay of R and their mass is larger than the "mass" of R, i.e. m, the LSP production 
could be sufficiently suppressed to make a reasonable contribution to dark matter. 

These and other manifestations of the considered scenario will be discussed elsewhere. 
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